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A new analytically and numerically manageable model collision operator is developed specifically 
for turbulence simulations. The like-particle collision operator includes both pitch-angle scattering 
and energy diffusion and satisfies the physical constraints required for collision operators; it con- 
serves particles, momentum and energy, obeys Boltzmann's _ff-theorem (collisions cannot decrease 
entropy), vanishes on a Maxwellian, and efficiently dissipates small-scale structure in the velocity 
space. The process of transforming this collision operator into the gyroaveraged form for use in 
gyrokinetic simulations is detailed. The gyroaveraged model operator is shown to have more suit- 
able behavior at small scales in phase space than previously suggested models. Model operators for 
electron-ion and ion-electron collisions are also presented. 



I. INTRODUCTION 

It has long been known that in many turbulent sys- 
tems the differences between vanishingly small dissipa- 
tion and neglecting dissipation completely are striking, 
and that this can be linked theoretically to the non- 
intcrchangeability of limits t ^ 00 and v 0, where 1^ 
is, e.g., viscosity, resistivity or collision frequency. Phys- 
ically, the dissipation is important in turbulence for the 
following reason. The fundamental property of turbu- 
lence is to transfer energy from scales at which it is in- 
jected into the system to scales where it is dissipated, 
leading to heating. When the dissipation coefficients are 
small, the system has to generate very fine-scale fluctu- 
ations in order to transfer the energy to scales at which 
dissipation becomes efficient. 

Because of Boltzmann's i?-theorem,— dissipation 
(meaning any effect that leads to irreversible heating) 
in kinetic plasmas is ultimately collisional, so the trans- 
fer of energy generally occurs in phase space — i.e., both 
in position and velocity space (see extended discussion of 
energy cascade in plasma turbulence in Ref. [H and ref- 
erences therein). There are a number of specific mecha- 
nisms, both linear and nonlinear, that give rise to phase- 
space mixing j^i^i^i^i^iii^ It is the resulting large gradients 
in the velocity space that eventually bring collisions into 
play however small the collision frequency is (such small- 
scale velocity-space structure has, e.g., been found and 
explicitly measured in gyrokinetic simulat ions^i^ i ^ ° 1 ^ ^ ) . 
Thus, in any plasma turbulence simulation, some effec- 
tive collisionality has to be present in order to smooth 
the small-scale structure in velocity space. 

Besides velocity-space smoothing, there is another 
key reason why collisions must be included. Colli- 
sions, through the dissipation of small-scale fluctuations 
in phase space, provide the physical link between ir- 
reversible plasma heating (macroscopic transport) and 



turbulence which enables the system to converge to a 
statistically steady state. Although it is possible for a 
collisionless simulation to temporarily achieve a quasi- 
steady state in macroscopic quantities, achieving a true 
steady state in the long-time limit requires some form of 
dissipation;^ While many simulations in plasma physics 
and neutral fluid dynamics have used numerical dissipa- 
tion, such as simple hyperdiffusion (or more sophisticated 
subgrid turbulence models in Large Eddy Simulations), 
to provide the dissipation needed for steady state, it is 
important to also be able to carry out direct numerical 
simulations, where the physical dissipation processes are 
explicitly resolved. This provides a valuable cross-check 
on simulations with numerical dissipation, and is useful 
as a standard by which one could search for optimal sub- 
grid models. 

Let us explain in more detail why collisions are im- 
portant for achieving the steady state. Consider the 
"(5/ kinetics," i.e., assume that it is physically reason- 
able to split the distribution function into a slowly (both 
spatially and temporally) varying equilibrium part and 
a rapidly varying fluctuating part: f = Fq + 5f. We 
further assume that Fq is a Maxwellian distribution, 
Fq = {no/n^^^v^^)exp{—v'^/v^Y^), where no is density, 
Vth = (2X0/"!)^/^ the thermal speed, Tq temperature and 
TO the particle mass. This will be the case if collisions are 
not extremely weak (for the weakly collisional formula- 
tion of Sf gyrokinetics, see Ref. [l3l ). One can show that 
the fundamental energy balance governing the evolution 
of the turbulent fluctuations is^i^iii^iHiiiii^ii^ 

where s is the species index, 5S = — JJ drdv 5f^ /2Fo is 



(1) 



2 



the entropy of the fluctuations, U = J dr {E"^ + B^)/87r 
is the energy of the (fluctuating) electromagnetic field, 
P is the input power (energy source of the turbulence), 
and C[(5/] is the linearized collision operator. In many 
types of plasma turbulence studied in fusion contexts, 
the input power P is proportional to the heat flux and 
it is the parameter dependence of the mean value of 
this quantity in the statistically stationary state that 
is sought as the principal outcome of the simulations. 
We can see immediately from the above equation that 
collisions (or some form of dissipation) are required to 
achieve such a steady state (as has been shown in numer- 
ical simulations^iiii^iii) and that in this steady state, P 
must be balanced on the average by the dissipation term. 

A key property of the collision operator required for 
this transfer of energy from turbulence to the equilibrium 
distribution to work correctly and, therefore, for the heat 
fluxes to converge to correct steady-state values, is that 
the collision term in Eq. ([1]) must be negative-definite: 

^['^■^l'^^^'" - (2) 

This ensures that heating is irreversible and that colli- 
sions cannot decrease entropy, the latter being the state- 
ment of Boltzmann's i/-theoremji While the heat fluxes 
might not be sensitive to the exact form of a model col- 
lision operator (within some range of models) at low col- 
lision frequency, any spurious sink of entropy may ad- 
versely affect the balance between turbulent fluxes and 
dissipation. Therefore it is clearly preferable that a 
model collision operator respect the i/-theorem, which 
has important physical consequences. Preserving the H- 
theorem may be even more important at higher collision 
frequencies. 

In view of the above discussion, we can formulate a rea- 
sonably restrictive set of criteria for any model collision 
operator: providing dissipation at small scales, obeying 
the i/-theorem [Eq. and also, obviously, conserv- 

ing particle number, momentum, energy, and vanishing 
on a (local; perturbed) Maxwellian distribution. Whilst 
these properties are analytically convenient, for numer- 
ical simulations the operator should also be efficiently 
implementablc and carry these properties (at least ap- 
proximately) over to the numerical scheme. 

The effect of small angle Coulomb collisions on an ar- 
birtrary distribution function was originally calculated 
by LandauF^^ In the Sf kinetics we would naturally con- 
sider the linearized Landau operatorJ^ However, it is 
sufficiently complex that a direct numerical convolution 
evaluation of it would exceed the limits on numerical re- 
sources that can be realistically expended on modelling 
the collisional physics of predominantly coUisionless plas- 
mas. Consequently several simplified model collision op- 
erators have been developed, both for analytical and com- 
putational convenience, that try to capture the qualita- 
tive essence, if not the quantitative detail, of the physics 
involved ]^°i^^'^^ This course of action is, indeed, emi- 
nently sensible: from Eq. ([T|), it seems plausible that. 



at least as far as calculating integral characteristics such 
as the turbulent fluxes is concerned, neither the exact 
functional form of the collision operator (provided it sat- 
isfies the criteria discussed above) nor the exact value of 
the collision frequency (provided it is sufficiently small) 
should be important. All we need is a physically reason- 
able dissipation mechanism. 

For these purposes, it has often been deemed suffi- 
cient to use the pitch-angle-scattering (Lorentz) opera- 
tor, sometimes adjusted for momentum conservation! 
However, in kinetic turbulence, there is no reason that 
small-scale velocity-space structure should be restricted 
to pitch angles. In fact, standard phase-mixing mecha- 
nisms applied to gyrokinetics produce structure in fy^'^ 
and there is also a nonlinear gyrokinetic phase mixing 
that gives rise to structure in v±, which may be an even 
faster and more efficient process.— Thus, a priori one 
expects to see small scales both in the pitch angle and in 
the energy variables (^ and v). It has, indeed, been con- 
firmed in simulationai^ that with only Lorentz scattering 
structure rapidly forms at the grid scale in energy. Thus 
a numerically suitable model collision operator should in- 
clude energy diffusioni^ 

In this paper, we propose such an operator (other op- 
erators including energy diffusion have been previously 
suggested j^^'^^ we include a detailed comparison of our 
operator with these in Appendix [Cjl . Our model opera- 
tor for like-particle collisions, including both pitch- angle 
scattering and energy diffusion and satisfying all of the 
physical constraints discussed above, is given in Sec. |TT] 
(the proof of the i7-thcorcm for it is presented in Ap- 
pendix In Sec. mil it is converted (gyroaveraged) 
into a form suitable for use in gyrokinetic simulations — 
a procedure that produces some nontrivial modifications. 
In Sec. IIVI we explain how interspecies collisions can be 
modelled in gyrokinetic simulations to ensure that such 
effects as resistivity are correctly captured. Section |V] 
contains a short summary and a discussion of the conse- 
quences of the work presented here. 

The anlytical developments presented in this paper 
form the basis for the numerical implementation of colli- 
sions in the publicly available gyrokinetic code GS2. This 
numerical implementation, as well as a suite of numeri- 
cal tests are presented in the companion paper, Ref. [lol 
(henceforth Paper II). 



II. A NEW MODEL COLLISION OPERATOR 

In this section, we present a new model collision op- 
erator for like-particle collisions that satisfies the criteria 
stated above. The interpecies collisions will be consid- 
ered in Sec. IIVI 

Let us start by introducing some standard notation. 
In discussing collision operators on phase space, we shall 
denote by r the position variable in physical space and 
use (u,^,!?) coordinates in velocity space, where v — \v\ 
is the energy variable, ^ = is the pitch-angle vari- 
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able, and it} the gyroangle about the equihbrium magnetic 
field. One can easily adapt the operators presented here 
to unmagnetized plasmas, but as we arc interested in gy- 
rokinetic plasmas, we shall concentrate on the strongly 
magnetized case. Taking the notation of Ref. [l^ as the 
standard, we introduce the normalized velocity variable 
X = v/vth and a set of velocity-dependent collision fre- 
quencies for like-particle collisions: 
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where = (2/^/??) e ^^dy is the error function, 

G{x) = [^{x) — x^'{x)]/2x^ is the Chandrasekhar func- 
tion, and v = -\/27rnog^ In ATq ^^^m~^/^ is the dimen- 
sional like-particle collision frequency (here In A is the 
Coulomb logarithm and q is the particle charge). Note 
that the two differential identities given in Eqs. ([5]) and 
([7]) will prove very useful in what follows. 

If one wishes to construct a model linearized collision 
operator, the following general form constitutes a natural 
starting point 



C[Sf] 



d_ 
dv 



D(v) 



d_5f 

dv Fa 



P[Sf]iv)Fo, (8) 



where the first term is the "test-particle" collision oper- 
ator and the second term the "field-particle" operator. 
Most model operators can be obtained by picking a suit- 
ably simple form for the velocity-space diffusion tensor 
D and the functional P, subject to the constraints that 
one chooses to impose on the model operator. 

In constructing our model operator, we retain the exact 
form of D for the linearized Landau collision operatorj^^ 



c[Sf] = i^oHSf]- 



)_d_ 

dv 



1 4 ^ d Sf 

2 " dv Fq 



-P[Sf]{v)Fo, 

where we have explicitly separated the energy-diffusion 
part (the second term) and the angular part (the first 



term), which includes pitch-angle scattering and is de- 
scribed by the Lorentz operator: 



m] = \ 



d ,2.96 f 1 dHf 
dC ^ ' d^ l-e dd'^ 



(10) 



Our modelling choice is to pick P to be of the form 



P[5f]{v) = + ^ Q[5f]. 



(11) 



One can view this prescription as first expanding P in 
spherical harmonics (one can easily show that they are 
eigenfunctions of the full field-particle operator) , reatain- 
ing only the first two terms, and then arbitrarily factor- 
izing the explicit v and 5f dependence of each harmonic. 
The functionals U[5f] and Q[5f] are mandated to have 
no explicit velocity dependence. In this ansatz the v de- 
pendence is chosen so that the final operator is self ad- 
joint and also to ensure automatic particle conservation 
by the field-particle operator: J P[5f]{v)FQdv = 0. In- 
deed the first term in Eq. (jlip gives a vanishing contribu- 
tion to this integral because it is proportional to v, and 
so does the second term because of the differential iden- 
tity given in Eq. The functionals U[Sf] and Q[Sf] 
are now uniquely chosen so as to ensure that the model 
operator conserves momentum and energy: a straightfor- 
ward calculation gives 



U[Sf] = K.vSfdv 1 1 {v/v,,,fv,Fodv, (12) 
Q[5f] = / v^i^ESfdv / f v^v/vti^fi^EFodv. (13) 



These are in fact just the standard correction ex- 
pressions used for the model pitch-angle-scattering 
operator-i2i^ and for more complex operators including 
energy diffusion.— 

Note that the numerical implementation of our colli- 
sion operator documented in Paper IliS achieves exact 
satisfaction of the conservation laws by choosing the dis- 
cretization scheme that exactly captures the differential 
identities Eqs. (O and ([7]) and the double integration by 
parts needed in deriving Eqs. ^2]) and p3)) . 

To summarize, we now have the following model oper- 
ator for like particle collisions: 



1 a /I 4 d 6f 

v^d~v i 2"" '^'i 



2v ■ U[6f] 



Fn 



v^ 



Q[6f]Fo, (14) 



th 
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where the functionals U[5f] and Q[Sf] arc given by Eqs. 



and (fT51) . The modelling choice of the field-particle 
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operator that we have made [Eq. (fTTj) ] means that, in 
order to compute our coUision operator, we have only 
to calculate definite integrals over the entirety of the ve- 
locity space — a significant simplification in terms of 
computational complexity and ease of use in numerical 
simulations (compared to computing convolutions over 
velocity space, see Paper IliS). 

As we have shown above, our operator conserves par- 
ticles, momentum and energy by construction. It is 
also not hard to see that it vanishes precisely when 
SJ/Fq = \,v,v^ and linear combinations thereof, i.e. if 
5f is a perturbed Maxwellian. From this and the fact 
that the operator is self adjoint, it can be easily shown 
that the operator only conserves particles, momentum 
and energy and that no spurious conservation laws have 
been introduced by our model. Because the model we 
have chosen retains the exact Landau test-particle op- 
erator , it provides velocity-space diffusion both in en- 
ergy and in pitch angle and will thus efficiently dissipate 
small-scale structure in velocity space. Finally, our model 
satisfies the iJ-theorem — this is proved in Appendix [XI 

This operator thus fulfills the criteria set forth in Sec. 
H] to be satisfied by any physically reasonable model op- 
erator. We now proceed to convert this operator into a 
form suitable for use in gyrokinetics. 



III. COLLISIONS IN GYROKINETICS 

The gyrokinetic theory is traditionally derived for a 
collisionless plasmaj^i^ However, as we have argued in 
Sec. m even when the collision frequency is small, colli- 
sions must be included in order to regularize the phase 
space and to ensure convergence of fluxes to statisti- 
cally stationary values. Mathematically, collisions can 
be included in gyrokinetics if the collision frequency is 
formally ordered to be comparable to the fluctuation 
frequency^ v ^ lo ^ fc||Wtii — the weakly collisional 
limit (collisionality larger than this leads simply to fluid 
equations). In practice, the collision frequency tends to 
be smaller than typical fluctuation frequencies, but this 
need not upset the formal ordering as long as it is not 
too small: the cases v ^ uj and u u can be treated as 
subsidiary limits 

Under the formal ordering v ^ uj ,ii is possible to show 
that the equilibrium distribution function (lowest order 
in the gyrokinetic expansion) is a MaxwellianH and the 
full distribution function can be represented as 



where Fq is a Maxwellian, ip the electrostatic potential (a 
fluctuating quantity) and h the (perturbed) distribution 
function of the paticle guiding centers. Here e — 
is the particle energy, ji = •mv\/2BQ the flrst adiabatic 
invariant, Bq the strength of the equilibrium magnetic 
field, R = r — p = r — bxv/Q. the guiding center po- 
sition, Vl the cyclotron frequency, and h = Bq/Bq. The 
gyrokinetic equation, written in general geometry and 
including the collision operator is then 



dh f dh c , , , , . 



dFo d{x)n 
de dt 



Bt 



{^o,(x>r} + CgkN,(16) 



where x = (y9 — • A/c the gyrokinetic potential, {x)r = 
(l/27r) J x{R + p) di} is an average over gyroangles hold- 
ing R fixed (the "gyroaverage" ) , vd is the guiding center 
drift velocity. 

The gyrokinetic collision operator CGK[/i] is the gy- 
roaverage of the linearized collision operator. The latter 
acts on the perturbed distribution h holding the parti- 
cle position r (not the guiding center R\) fixed. This 
nuance must be kept in mind whilst working out the ex- 
plicit form of CQK[h{R)] from the unaveraged linearized 
operator C[h{r — p)]. 

Let us restrict our consideration to local simulations, 
which are carried out in a flux tube of long parallel extent, 
but short perpendicular extent. In such simulations, one 
assumes that the equilibrium profiles are constant across 
the tube, but have non-zero gradients across the tube 
so as to keep all the appropriate drifts and instabilities. 
This permits one to use periodic boundary conditions 
and perform the simulations spectrally perpendicular to 
field lines Thus 



(17) 



/ 



h{t, R, ^, e) 



(15) 



where I is a coordinate along the field line and the Fourier 
transform is understood to be only with respect to the 
perpendicular components of R, i.e., k = k±. Treat- 
ing the perpendicular coordinates spectrally confines all 
dependence on the gyroangle to the exponent, thus 
we can compute the gyroangle dependence explicitly and 
carry out the gyroaveraging of the collision operator in a 
particularly transparent analytical way)^i^ 



CgkN = ( C 



ik-R 



hk 



R 



k 



(18) 
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where p = b y. v^/Vl. Thus, in Fourier space 

CGK[hk] = {e'''-PC[e-'''-Phk\), (19) 

where (...) refers to the expUcit averaging over the 
dependence. Some general properties of this operator 
are discussed in Appendix B of Ref. [Sj. 

We now apply the general gyroaveraging formula Eq. 
(fT9ll to our model operator given by Eq. (fT4|) . The gyroki- 
netic transformation of variables (r , w, ^, ■&) — > (/?, ^, e, d) 
mixes position and velocity space. However, in the colli- 



sion operator, to the lowest order in the gyrokinetic ex- 
pansion, we can neglect the spatial dependence of ^ that 
comes via the equilibrium magnetic field B^ir) and thus 
use the (w,^) velocity variables. After some straightfor- 
ward algebra, which involves converting velocity deriva- 
tives at constant r to those at constant R and evaluating 
the arising gyroaverages as detailed in Appendix [B] we 
arrive at the following model gyrokinetic collision opera- 
tor 



where p = Vth/^ is the thermal Larmor radius (not to be confused with the velocity-dependent p), a = kj_vj_/D,, Jq 
and Ji are Bessel functions and 



U±[hk] = 2] ^sV^Ji{.a)hkdv j J (v/vth) K-iFadv, (21) 
U\\[hk] = ^ j VsV\\Jo{a)hkdv j J {v/vthf i^sFq dv, (22) 
Q[hk] = I v^VEJG{a)hkdv v'^ (v/vthf veFq dv. (23) 



Note that since the position and velocity space are mixed 
by the gyrokinetic transformation of variables, R = r—p, 
the collision operator now contains not just pitch-angle 
and V derivatives but also a spatial perpendicular "gy- 
rodiffusion" term. 

It important to make sure that the operator we have 
derived behaves in a physically sensible way in the long- 
and short-wavelength limits. When k±p ^ 1, all the 
finite-Larmor-radius (FLR) effects, including the gyrod- 
iffusion, disappear and we end up with pitch-angle scat- 
tering and energy diffusion corrected for energy and par- 
allel momentum conservation — the drift-kinetic limit. 
In the opposite limit, k±p ^ 1, we can estimate the 
behavior of our operator by adopting the scaling of the 
velocity derivatives based on the nonlinear perpendicu- 
lar phase mixing mechanism for gyrokinetic turbulence 
proposed in Ref. 0: this produces velocity-space struc- 
ture with characteristic gradients vthd/dv± ^ k±p (see 
also Refs. HQ). With this estimate, we see that all the 
field-particle terms in the operator are subdominant by 
a factor of {k±p)~^. Thus the operator reduces to the 
gyrokinetic form of the test-particle Landau operator in 
this limit. All diffusive terms are also equally large in this 
scaling, supporting our supposition that energy diffusion 
needs to be included. These considerations give us some 



confidence that we correctly model the diffusive aspects 
of the coUisional physics in a short-wavelength turbulent 
regime. Indeed, if one applies the same estimates to the 
full linearized Landau operator, the Rosenbluth poten- 
tials of the perturbation arc small when k±p ^ 1 be- 
cause they are integrals of a rapidly oscillating function, 
so the dominant effect does, indeed, come entirely from 
the test-particle part of the operator. 

The gyrokinetic if-theorem, which has to be satisfied 
in order for heating and transport to be correctly calcu- 
lated, is given by^ii^ 

JI ^^CGK[h]dRdv<Q. (24) 

The gyrokinetic collision operator given by Eq. (|20p re- 
spects this inequality, as can either be shown directly 
from Eq. (pO|) (analogously to the proof in Appendix \^ 
or inferred from Eq. ^ by transforming to gyrokinetic 
variables. The operator also manifestly diffuses small- 
scale structure both in velocity and in (perpendicular) 
position space. 

How to express the conservation-law tests upon this 
operator is a somewhat subtler question. This is because 
after the gyroaveraging has been done, one cannot ex- 
plicitly separate the position- and velocity-space dynam- 
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ics in the gyrokinetic phase space. However, it is still 
possible to express the evolution of particles, momentum 
and energy as local conservation laws. Let us take the 
velocity moments of the gyrokinetic equation Eq. (jl6p 
corresponding to these conserved quantities. These are 
evaluated at constant position r = R+p. Defining (. . 
as the average over gyroangles while holding r fixed and 
using Eq. (|18p , we arrive at the following evolution equa- 
tion for the conserved moments: 




-ik p \ /Jk p 



v\\b I {C[hk])dv - ik-T'^ 



(coll) 



(25) 



where f '^''^ denotes the fluxes arising from the terms in 
Eq. (fT6)) other than the collision operator. The mo- 
ments of the gyroaveraged collision operator have been 
expressed as the moments of the same operator formally 
taken at fc • p = (i.e., dropping all FLR contribu- 
tions: {C[hk]) instead of (e*'=PC[e-''=P/ifc])) plus the 
divergence of the coUisional flux arising from the finite- 
Larmor-radius part of Cgk[/j]- This representation was 
achieved by expanding all gyrophase factors in Eq. ([25]) 
into infinite Taylor series and noticing that they take the 
forme±''='' = l±ik- p~ {l/2){k- pf+ ■ ■ ■ = •[...], 
where the square brackets contain the rest of the series. 
Without the FLR terms, the particle, momentum and 
energy moments of the gyrokinetic collision operator van- 
ish, 




{C[hk]) dv = 0, 



(26) 



and Eq. (|25p represents the local conservation law for 
these quantities, with the fluxes containing both coUi- 
sionless and coUisional contributions: T = T^^^i -\-T('^°^^'> ^ 
Thus, a conservative numerical implementation of the gy- 
rokinetic collision operator can be achieved if Eq. ([^5]) is 
hard-wired into the numerical sheme.— How to do this is 
explained in Paper H^iS where we also demonstrate the 
correct performance of our model operator on a number 
of test problems. 



IV. INTERSPECIES COLLISIONS 

Let us now turn to the collisions between different 
species and focus on a plasma containing only electrons 
and one species of ions with a mass ratio m^/rrii <IC 1. 
The smallness of the mass ratio allows for a significant 
simplification of the interspecies collision terms. Since 
ion-electron collisions are subdominant to the ion-ion 
oneSfi^ Vie/vii ~ {rrie/miy^^ , to lowest order in the mass 
ratio we can neglect the ion-electron collisions and the ion 
collisions can be modelled using the like-particle opera- 
tor proposed above [Eq. (HOI)]- The situation is different 
for the electron- ion collisions, which are same order in 
mass ratio as the electron-electron collisions)^ Pei ~ ^^ee- 
Thus, the full electron collision operator has two parts: 



C[dfe]^Cee[S.fe]+Ce^dSf]. 



(27) 



The electron-electron operator Cee[Sfe] can be modelled 
by the like-particle operator proposed above [Eq. (dH)], 
the electron-ion collision operator can be expanded in the 
mass ratio and to two leading orders reads^^ 



2?? • 11 



''the 



(28) 
(29) 



where i/^i = v^TTTioi-Z^^e^lnATQg'^^^me is the dimen- 
sional electron- ion collision frequency, Z = qi/e, e is the 
electron charge, L is the Lorentz operator given by Eq. 
(fTU)l . and 



— I vSfi dv 



(30) 



is the ion flow velocity. Thus, the electron- ion collisions 
are correctly modelled to lowest order in the mass ra- 
tio by electron pitch-angle scattering off static ions, with 
electron drag against the bulk ion flow as a first order 
correction to this. 

Note that the drag term is necessary to correctly cap- 
ture electron-ion friction and hence resistivity; failure to 
include it leads to incorrect results, with mean electron 
momentum relaxed towards zero rather than towards 
equality with the mean ion momentum. The need to in- 
clude this first-order effect stems from the fact that whilst 
the effect this has on the mass flow is small (mass is pre- 
dominantly carried by the ions) there is a large effect on 
the current. Including the small correction to the mean 
ion motion due to friction on electrons results in a small 
correction to both the current and the mass flow. How- 
ever, these first-order effects such as the small slow col- 
lisional change in the mean ion momentum, are formally 
of the same order as the drag terms in the electron-ion 
collision operator. We will, therefore, keep the first-order 
correction to the ion collision operator. 

Taking the lowest-order contribution to the linearized 



7 



ion-clcctron operator— 

Qe[5/]--^-^, (31) 

where Rie is the ion-electron friction force. Since Rie = 
—Ret = — J iTT-ef^CeilS fe\dv , the ion-clcctron collision op- 
erator is expressed in terms of the perturbed electron dis- 
tribution function: using lowest-order term in Eq. (j29p 
(the pitch-angle scattering), we find, to lowest order in 
the mass-ratio expansion, 

C^e[Sfe] = ■ [ v' m,vf,{v')5 fe{v')dv' . (32) 

It is now easy to see that formally the ion-electron col- 
lisions must be kept in order for the i?-theorem to be 
satisfied. Indeed, for the interspecies collisions, the H- 



theorem is written as follows 

/ 1 ^^-t'^/]^^'^^ + // ^^^C,A5f]dvdr < 

(33) 

(i.e., the interspecies terms in Eq. ^ are nonpositive) . 
We can see immediately that the pitch-angle-scattering 
part of the electron-ion operator automaticlaly satisfies 
this, while the contribution from the drag term in Eq. 

is exactly cancelled by the contribution from the 
ion-electron operator given by Eq. (|32p . 

Let us now work out the gyrokinetic interspecies colli- 
sion operators. Performing the conversion of the electron- 
ion operator [Eq. ((29l) ] to the gyroaveraged form in a way 
entirely analogous to what was done in Sec. IIIII and Ap- 
pendix |B1 we get 



Cgk [^efej 



"1 d_ 

2 d£, 



(1-e 



2\ dhek 



9( 



ZrUe v]_ Ji(ae) 



2w||Jo(ae)u|| 



ik 



Oe 



"the 



2 „2 



Foeklp: 



""the 



2v'^^J,{a[) 



hik{v')dv' 



"tht 



(34) 



where 

u\uk = — [ vuJo{ai)hikdv (35) 
noi J 

and fls = k±v±/fls for species s and the rest of the nota- 
tion as in previous sections, except with species indices 
reintroduced. 

Let us estimate the size of the four terms in Eq. ([34]) 
at the ion (long) and electron (short) scales. The first 
term (pitch-angle scattering) is always important. At 
the ion scales, k±pi ^ 1, the third term (parallel ion 
drag) is equally important, while the second term (elec- 
tron gyrodiffusion) and the fourth term are subdom- 
inant by a factor of me/rrii. At the electron scales, 
k±Pe ~ 1, the pitch-angle scattering and the electron 
gyrodiffusion (the first two terms) are both important. 



Since at these scales kj_pi ^ [mi/raeY/'^ ^ 1, the third 
and fourth terms are subdominant by a factor (resulting 
from the Bessel functions under the velocity integrals) of 
l/V^'J-Pi ~ (fTie/'Tii)^/'*. In fact, they are smaller than 
this estimate because at these short wavelengths, the ion 
distribution function has small-scale structure in veloc- 
ity space with characteristics scales <5w_L/fthi ^ ^/k±pi, 
with leads to the reduction of the velocity integrals by 
another factor of 'Vj \fkZpi- Thus, at the electron scales, 
the third and fourth terms in Eq. ([34]) are subdominant 
by a factor of (mg/mi)^/^. 

These considerations mean that the fourth term in 
Eq. p4|) is always negligible and can safely be dropped. 
The full model gyrokinetic electron collision operator is, 
therefore. 



1 d_ ^2-. ^Kk 

2 dC ' 



^ ^the 



V 
,,2 



k\p\hf,k 



2t;|| Jo(ae)M||j 



Oe 



(36) 



where the electron-electron model operator CQY}^ek\ is a form that does not contain an explicit dependence on 
given by Eq. ([20]) and u\\ik by Eq. (|35|) . Note that since 
enoe(u||i — U||e) = 3\\ is the parallel current, the parallel 
Ampere's law can be used to express in Eq. (|36p in 
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the ion distribution function: 

If c 

uuk = / v\\jQ{ae)hek dv + 



tation of the electron operator, detailed in Paper II 
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fciAiifc. (37) 

Performing a completely similar calculation for the ion- 
electron operator [Eq. ([5^ ]. we obtain the full model 



This turns out to be useful in the numerical implemen- gyroaveraged ion collision operator: 
I 



CoKlhik] = CcKlf^ik] — 



F< 



Oi 



Oi 



wii Jo(ai) / v[jQ{a'^)mev''j^{v')hek{v')dv' , 



r 



(38) 



where the ion-ion model operator CQy^lh^f,] is given by 
Eq. (EOl). 



V. SUMMARY 

Thus in Sec. H] we have argued the necessity of dissipa- 
tion in turbulence simulations, justified the direct mod- 
elling of collisions in order to provide such dissipation and 
postulated a set of constraints for a physically reasonable 
model collision operator. Previously used model opera- 
tors mostly do not contain energy diffusion and were thus 
deemed unsatisfactory for these purposes. Of the exist- 
ing model operators that contain energy diffusion, two 
are detailed in Ref. [13 and Ref . [2l|. The former however 
does not satisfy the 7?-theorem [Eq. ^] and the latter 
incorrectly captures the smallest scales. These problems 
are demonstrated and discussed in detail in Appendix ICl 

In Sec. HI] we presented a new operator [Eq. (fT4|) ] that 
successfully introduces energy diffusion whilst maintain- 
ing the i7-Theorem and conservation laws, thus satis- 
fying the conditions set forth in the introduction. This 
operator is then transformed into gyrokinetic form in Sec. 
mil correctly accounting for the gyrodiffusive terms and 
FLR effects [Eq. ((20)) ]. In order to provide a complete 
recipe for modelling the coUisional effects in simulations, 
the same gyroaveraging procedure is applied in Sec. lIVI to 
electron-ion and ion-electron collisions, somewhat sim- 
plified by the mass-ratio expansion [Eq. ([36]) and Eq. 
([55]) ]. This leaves us with a complete picture of colli- 
sions in gyrokinetic simulations, capturing gyrodiffusion, 
resitivity and small-scale energy diffusion. 

When we discussed the gyroavergaing procedure in 
Sec. ini] we presented the specific case of the applica- 
tion to Eulcrian flux-tube 5f gyrokinetic simulationsi^i^ 
However, the form presented in Eq. p4p is suitable for 
inclusion in most 5f kinetic systems and even amenable 
to use in Lagrangian codes by applying the methods of 
Rcfs. |28i ori2^ to the gyroaveraged operator given by Eq. 
(|20)) . Indeed, by suitable discretization of the gyroav- 
eraging procedure^! it would also be usable in a global 
Eulerian code. 

We conclude by noting that the final arbiter of the 
practicality and effectiveness of this collision model is 



the numerical implementation and testing performed in 
Paper II,— where our operator is integrated into the GS2 
code. The battery of tests shows that our operator not 
only reproduces the correct physics in the weakly coUi- 
sional regime but even allows a gyrokinetic code to cap- 
ture correctly the coUisional (reduced-MHD) limit. 
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APPENDIX A: PROOF OF THE //-THEOREM 
FOR EQ. ([14]) 

In the case of the expansion / = Fo + (5/ about a 
Maxwellian the entropy generation by like particle colli- 
sions takes the form 



dS_ 



d_ 

dt 



f In / dv dr 
fC[fFo]dvdr>Q, (Al) 



where we use the compact notation / = 5f/Fo. The 
statement of the //-theorem is that the right-hand side 
of Eq. (|Aip is nonnegative and that it is exactly zero 
when Sf is a perturbed Maxwellian. 

We represent / as a Cartesian tensor expansion (or 
equivalently spherical harmonic expansion) in velocity 
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space: 

fir,v)^ foir,v)+v fi{r,v)+R[f]ir,v), (A2) 

where R[f] comprises the higher order terms. It is then 
possible to recast the statement of the ff- Theorem in 
terms of this expansion using linearity of the model col- 
lision operator C [Eq. (fT4ll ] . orthogonality of the expan- 
sion and the fact that spherical harmonics are eigenfunc- 
tions of the Lorentz operator L. By construction, R[f] 
satisfies / R[f]Fodv — and J vR[f]Fodv = 0, from 
which it follows that R[f] docs not contribute to the field- 
particle parts of the model operator: U[R[f]Fo] = and 
Q[R[f]Fa] = 0. Substituting Eq. jJS]) into the right- 
hand side of Eq. (jAip . where the operator C is given by 
Eq. p4p . and integrating by parts those terms involving 
derivatives of we find that they all give nonnegative 
contributions, so we have 



- / fC[fFo]dv>ao + ai, 



where 



(A3) 

(A4) 
(A5) 



In order to prove the _ff-thcorcm, it is now sufficient to 
show that (To > and cti > 0. 

Starting with ctq and using Eq. (HU, we integrate over 
angles and use the differential identity given in Eq. ([7]) 
to express the term containing Q: 



c^o = - y foC[foFo]dv, 
ai = - f V- fiC[v- AFo]dv. 



CTo = -27r / /o 



d_ 

dv 



dv V i4 



dv. 
"(A6) 

Using the aforementioned identity again in the expression 
for Q [Eq. (|13p ] and integrating by parts where oppor- 
tune, we get 



ctq = 47r 




5 r J 
-—V lynFodv 

ov " 




v^VEF{)dv 



(A7) 



It is easy to see from the Cauchy-Schwarz inequality that 



^v^y^^Fodv] < 




v'^h'^Fodv I v^v^Fodv. 

(A8) 



Using this in the second term of Eq. (jA7[) . we infer 

2 

CTo > 47r 




v'^i^^^Fodv 

\ - Jv'^u^^Fodv j Jv'^iyEFodv^ = 0, 



(A9) 



where to prove that the right-hand side vanishes, we 
again used the differential identity given in Eq. ([7]) and 
integrated by parts. Thus, we have proved that ctq > 0. 

Turning now to cti [Eq. (|A5p ]. using Eq. (|14p . and in- 
tegrating by parts where opportune, we get 



0-1 = I [v- fifvoF^dv I { • fi ] v^v\\Fodv 



3v, 



th 



XX ■ fi^sFodv 




x^i'sFodv, 



(AlO) 



where we have used the standard notation that x = v /vth 
and X = v/vth- Integrating over angles and using the 
simple identity a ■ J vvdil = (47r/3)a, where i) = v/v 
and a is an arbitrary vector, we have 



CTl 



47rut\ 



1 



Ifil^x^i^oFodx + - 



d , 

ox 



fix i^sFodx 




X UgF^dx 



x^^v^F^dx 



(All) 



Once again applying the Cauchy-Schwarz inequality, we 
find that 



fix*i^sFodx 



< 



fi 



x'^UgF^dx I x i^gFodx 



Using this in the last term in Eq. (jAlip . we get 
An 



(^1 > Y^'th 



/i 



x^Ai^Fodx 



d , 

ox 



(A12) 



(A13) 



X v/nF^dx 



where Av \s defined in Eq. Upon using the differ- 
ential identity given in Eq. ([6]) to express l^v in the first 
term of the above expression and integrating the result- 
ing expression by parts, we finally obtain 





dfi 




dx 



x^v\\Fodx > 0. 



(A14) 



We now consider when these inequalities becomes 
equalities, i.e., when the right-hand side of Eq. (|Aip 
is zero. Firstly, this requires dR[f]/d^ = and thus 
i?[/] = 0, so / contains no 2nd or higher-order spherical 
harmonics. Secondly, cto = if cither /o is independent 
of V or we have equality in the invocation of the Cauchy- 
Schwarz inequality [Eq. (jASP J. which occurs if /q oc u^. 
Similarly cti = iff /i is independent of v. Thus, the 
right-hand side of Eq. (|Aip vanishes iff / oc i.e., 
= JFo is a perturbed Maxwellian. 

This completes the proof of the _ff-theorem for our 
model operator. 
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APPENDIX B: GYRO AVERAGING 

To transform the derivatives in Eq. (|14p from the orig- 
inal phase-space coordinates {r,v,^,'&) to the new coor- 
dinates (i?, w, ^, we require the following formulae: 

where p = b x v±/^l. In Fourier-transformed perpen- 
dicular guiding center variables, we can replace in the 
above formulae {d/dR)^ ik, where k = k±. It is con- 
venient to align (without loss of generality) the 79 = 
axis with k, so we have v± ■ k = /cxi'v'T^^^cos'i? and 
p ■ k ^ —k±vy^l — ^2 sin-iJ. Using the above formulae, 
we gyroaveragc the Lorentz operator in Eq. (|14p : 

where we have used {pp)ji = {v±^Vi_)jj^ = (1/2)/. 
Note that both the terms containing ^ and d deriva- 
tives in the original operator Eq. (|10p produce non- 
zero gyrodiffusive contributions [the second term in 
Eq. (jB4[) ]. Another such gyrodiffusive term, equal 
to — [z;'^(l — ^^)/4r2] fcj^/ifc, arises from the energy- 
diffusion part of the test-particle operator in Eq. 
Collecting these terms together and defining the thermal 
Larmor radius p = wth/^^, we arrive at the gyrodiffusion 
term in Eq. pT]) . 



It remains to gyroaveragc the field-particle terms. For 
the energy-conservation term [Eq. (jl3p ] we have 

{e'^-PQle-^'^-Phu]) - (e"^-'') Qie-'^^ Phkl (B5) 

where 

J v^iyE{e-"'-'')hkdv ^ Jv^ {v/vthfvEFodv. 



Note that the integration in Q only affected e 
hence the above expression. Using the standard Bessel 

function identity^ / e*""'"''di9 27rJo(a), we find 



^g-ifc p^ — o/o(a), where a = k±v±/n. Substituting 
this into Eqs. (|B5p and (|B6p . we arrive at the energy- 
conservation term in Eq. (|20p . where the expression in 
the right-hand side of Eq. (|B6p is denoted Q[hk] [Eq. 

mi 

The momentum-conserving terms are handled in an 
analogous way: details can be found in Appendix B of 
Rcf.[g, where a simpler model operator was gyroaveraged. 
APPENDIX C: COMPARISON WITH 
PREVIOUS MODEL OPREATORS 

In order to compare and contrast with previously sug- 
gested operators that do include energy diffusion, we first 
rewrite in our notation the operator derived by Catto and 
Tsang [Eqs. (14) and (16) in Ref.lH], 



I 



This operator, whilst it conserves particle number, mo- 
mentum and energy, neither obeys the _ff-Theorem nor 
vanishes on a perturbed Maxwellian. 

The latter point can be demonstrated most easily by 
letting Sf ~ x^Fq, where x = w/wth- This 5f is pro- 
portional to a perturbed Maxwellian with non-zero 5n 
and 5T. We can then evaluate the test-particle and field- 
particle parts of the operator to find 



1 9 /I 4 d 5f\ Id,. . 2 
ov \2 " ov to J x'^ ox 

(C2) 



and 

00 — 
j x^i'sSfdv = J x^VEe~^ dx = y — noh'. (C3) 


Substituting into Eq. (jCip . we get 

CcT[Sf]^-x^i^EFo + ^(^x^-^^i^Fo, (C4) 

which is non-zero despite Sf being a perturbed 
Maxwellian. 

In order to show that the iJ-thcorem can be violated 
by the operator Eq. (jCip . let us consider a perturbed 
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distribution function of the form 5f = x^Fq. Then 
"3 , 



CcT[Sf] 



ve) 



Fo, (C5) 



so the entropy generation is, 

^^-JJy^C'^^ [-^/l dvdr = - A(32-21^/2) i^V < 0, 

(C6) 

where V is the volume of the system. The above expres- 
sion is negative, which breaks the iJ-theorem and pro- 



duces unphysical plasma cooling for the particular form 
of the perturbed distribution function that we have ex- 
amined. 



The second case we examine here is the sequence of 
operators derived by Hirshman and Sigmarj2i The gen- 
eral operator proposed by these authors is given in their 
Eq. (25). In our notation, we rewrite here the = 0, 
I = 0,1 restriction of the like-particle form of their oper- 
ator with Ai/ set to for simplicity (this docs not affect 
the discussion that follows): 



2 -TO, 



(C7) 



where U and Q arc defined by Eqs. p2p and (|13p . The primary concern here comes from the angle averaging operation 
in the energy-diffusion part of the operator. Firstly, the energy diffusion only acts on the spherically symmteric (in 
velocity space) part of the perturbed distribution function. However, there is no reason why there cannot arise 
perturbations that have very large energy derivatives but angle-average to zero (for example, 5f cx i^). Clearly, such 
perturbations will not damped correctly. Secondly, upon conversion to gyrokinetic coordinates and gyroaveraging (see 
Sec. mil and Appendix [B]) , the operator becomes 



HS,GK 



[hk] = J^d(w) 



ld_ 
2d^ 



dhk 



th 



v±Ji{a)U±_ [hk] + W|j Jo(a)J7|[ [h^ 



Joja) d 
dv 



1 4 p a 

2- -11^0^ 



2TTFn 



(C8) 



Ma)Q[hk]Fo, 



where the conservation functionals U±, U\\ and Q are 
the same as defined in Eqs. (PT|) - (|^5|) . The immcdiatly 
obvious problem is that the angle averaging has intro- 
duced two new Bessel functions into the energy diffusion 
term. The energy diffusion is therefore supressed by one 
power of k±p in the limit k^^p ^ 1, while it is precisely 
in this limit that we expect the small-scale structure in 
the velocity space to be particularly important i^i^ This 
means that the energy cutoff in phase space is artificially 
pushed to smaller scales and one might encounter all the 
problems associated with insufficient energy diffusioniiS 



While, for the reasons outlined above, we expect the 
Hirshman-Sigmar operator not to be a suitable model for 
collisions, we would like to note that for many purposes 
the Hirshman-Sigmar operators are superior to the model 
operator we presented in Sec. |TT1 Taken as a sequence, 
they provide a rigorous way of obtaining classical and 
neoclassical transport coefficients to any desired degree 
of accuracy, and it is relatively easy to solve the Spitzer 
problem for them, while the Spitzer functions for our 
operator are hard to find analytically. 
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